

* This file is a re-replication of Gibler, Miller, and Little replication of Braithwaite and Lemke (2011). 
* This is the second of two Braithwaite and Lemke (BL) replication files.
* It focuses on replicating Braithwaite and Lemke's statistical analysis.

clear
capture log close
cd "C:\Users\rum842\Dropbox\3.4 MID McManus\ISQ Response\Final Submission\Replication Files" //Insert your own directory here.
log using "Braithwaite and Lemke Analysis Replication.log", replace
set more off

* Merge in variables created in the previous do file.

use "B&L_CMPS.dta", clear
* Note: This is the BL replication dataset from Gibler's dataverse, but we added the variables ccodea and ccodeb and dropped some estimated variables.

* The dyadid is missing in some cases, so we fix this first.
gen ccodehigh=ccodea if ccodea>ccodeb
replace ccodehigh=ccodeb if ccodeb>ccodea
gen ccodelow=ccodea if ccodea<ccodeb
replace ccodelow=ccodeb if ccodeb<ccodea
replace dyadid=ccodelow*1000+ccodehigh

merge m:1 dyadid year using "BL Dyadic MIDs 3.0.dta"
drop if _merge==2
drop _merge
replace onset_mid3=0 if onset_mid3==.

merge m:1 dyadid year using "BL Dyadic MIDs 4.3.dta"
drop if _merge==2
drop _merge
replace onset_mid4=0 if onset_mid4==.

merge m:1 dyadid year using "BL Dyadic MIDs GML ISQ.dta"
drop if _merge==2
drop _merge
replace onset_gml1=0 if onset_gml1==.

merge m:1 dyadid year using "BL Dyadic MIDs GML 2.1.dta"
drop if _merge==2
drop _merge
replace onset_gml2=0 if onset_gml2==.

* Labeling
label var thom_rival "Rivalry"
label var min_min "Both Minors"
label var joint_democ "Joint Democracy"
label var jointsatis "Joint Satisfaction"
label var terrdisp "Territorial Claim"
label var powerratio "Power Preponderance"
label var defense "Allied" 
label var territory "Territorial Dispute"
label var contiguous "Contiguous"				


****************************************  ANALYSIS  **************************************************************************

* ANALYSIS WITH MID 3.0 DATA (ORIGINALLY USED BY BL)

* In addition to running the same regressions as BL, we add some matrix notation to help us make the figures.

* MID 3.0 Recip Model

heckprob recip_mid3 joint_democ territory jointsatis powerratio defense , ///
	sel(onset_mid3 = contiguous thom_rival min_min joint_democ  jointsatis ///
	terrdisp powerratio defense peaceyears pax2 pax3) cluster(dyadid)
eststo m1_bl
***
matrix m1_bl_out = J(5,2,.)
matrix coln m1_bl_out = coef se
matrix rown m1_bl_out = joint_democ territory jointsatis powerratio defense
matrix beta = e(b)
matrix cov = e(V)

foreach i of num 1/5 {
	matrix m1_bl_out[`i',1] = beta[1,`i']
	matrix m1_bl_out[`i',2] = sqrt(cov[`i',`i'])
}

matrix m1_bl_sel = J(9,2,.)
matrix coln m1_bl_sel = coef se
matrix rown m1_bl_sel = contiguous thom_rival min_min joint_democ jointsatis ///
						terrdisp powerratio defense Rho
*Pull each variable, plus Rho
foreach i of num 7/14 {
	matrix m1_bl_sel[`i'-6,1] = beta[1,`i']
	matrix m1_bl_sel[`i'-6,2] = sqrt(cov[`i',`i'])
}
matrix m1_bl_sel[9,1] = beta[1,19]
matrix m1_bl_sel[9,2] = sqrt(cov[19,19])
***

* MID 3.0 Use of Force Model

heckprob force_mid3 joint_democ territory  jointsatis powerratio defense , ///
	sel(onset_mid3 = contiguous thom_rival min_min joint_democ  jointsatis ///
	terrdisp powerratio defense peaceyears pax2 pax3) cluster(dyadid)
eststo m2_bl
***
matrix m2_bl_out = J(5,2,.)
matrix coln m2_bl_out = coef se
matrix rown m2_bl_out = joint_democ territory jointsatis powerratio defense
matrix beta = e(b)
matrix cov = e(V)

foreach i of num 1/5 {
	matrix m2_bl_out[`i',1] = beta[1,`i']
	matrix m2_bl_out[`i',2] = sqrt(cov[`i',`i'])
}

matrix m2_bl_sel = J(9,2,.)
matrix coln m2_bl_sel = coef se
matrix rown m2_bl_sel = contiguous thom_rival min_min joint_democ jointsatis ///
						terrdisp powerratio defense Rho
*Pull each variable, plus Rho
foreach i of num 7/14 {
	matrix m2_bl_sel[`i'-6,1] = beta[1,`i']
	matrix m2_bl_sel[`i'-6,2] = sqrt(cov[`i',`i'])
}
matrix m2_bl_sel[9,1] = beta[1,19]
matrix m2_bl_sel[9,2] = sqrt(cov[19,19])
***


* MID 3.0 Mutual Use of Force Model

heckprob mut_force_mid3 joint_democ territory  jointsatis powerratio defense  ///
	,sel(onset_mid3 = contiguous thom_rival min_min joint_democ  jointsatis ///
	terrdisp powerratio defense peaceyears pax2 pax3) cluster(dyadid)
eststo m3_bl
***
matrix m3_bl_out = J(5,2,.)
matrix coln m3_bl_out = coef se
matrix rown m3_bl_out = joint_democ territory jointsatis powerratio defense
matrix beta = e(b)
matrix cov = e(V)

foreach i of num 1/5 {
	matrix m3_bl_out[`i',1] = beta[1,`i']
	matrix m3_bl_out[`i',2] = sqrt(cov[`i',`i'])
}

matrix m3_bl_sel = J(9,2,.)
matrix coln m3_bl_sel = coef se
matrix rown m3_bl_sel = contiguous thom_rival min_min joint_democ jointsatis ///
						terrdisp powerratio defense Rho
*Pull each variable, plus Rho
foreach i of num 7/14 {
	matrix m3_bl_sel[`i'-6,1] = beta[1,`i']
	matrix m3_bl_sel[`i'-6,2] = sqrt(cov[`i',`i'])
}
matrix m3_bl_sel[9,1] = beta[1,19]
matrix m3_bl_sel[9,2] = sqrt(cov[19,19])
***


* MID 3.0 Fatal Model

heckprob fatal_mid3 joint_democ territory  jointsatis powerratio defense , ///
	sel(onset_mid3 = contiguous thom_rival min_min joint_democ  jointsatis ///
	terrdisp powerratio defense peaceyears pax2 pax3) cluster(dyadid)
eststo m4_bl
***
matrix m4_bl_out = J(5,2,.)
matrix coln m4_bl_out = coef se
matrix rown m4_bl_out = joint_democ territory jointsatis powerratio defense
matrix beta = e(b)
matrix cov = e(V)

foreach i of num 1/5 {
	matrix m4_bl_out[`i',1] = beta[1,`i']
	matrix m4_bl_out[`i',2] = sqrt(cov[`i',`i'])
}

matrix m4_bl_sel = J(9,2,.)
matrix coln m4_bl_sel = coef se
matrix rown m4_bl_sel = contiguous thom_rival min_min joint_democ jointsatis ///
						terrdisp powerratio defense Rho
*Pull each variable, plus Rho
foreach i of num 7/14 {
	matrix m4_bl_sel[`i'-6,1] = beta[1,`i']
	matrix m4_bl_sel[`i'-6,2] = sqrt(cov[`i',`i'])
}
matrix m4_bl_sel[9,1] = beta[1,19]
matrix m4_bl_sel[9,2] = sqrt(cov[19,19])
***


* MID 3.0 High Fatality Model

heckprob high_fatal_mid3 joint_democ territory  jointsatis powerratio defense , ///
	sel(onset_mid3 = contiguous thom_rival min_min joint_democ  jointsatis  ///
	terrdisp  powerratio defense peaceyears pax2 pax3) cluster(dyadid)
eststo m5_bl
***
matrix m5_bl_out = J(5,2,.)
matrix coln m5_bl_out = coef se
matrix rown m5_bl_out = joint_democ territory jointsatis powerratio defense
matrix beta = e(b)
matrix cov = e(V)

foreach i of num 1/5 {
	matrix m5_bl_out[`i',1] = beta[1,`i']
	matrix m5_bl_out[`i',2] = sqrt(cov[`i',`i'])
}

matrix m5_bl_sel = J(9,2,.)
matrix coln m5_bl_sel = coef se
matrix rown m5_bl_sel = contiguous thom_rival min_min joint_democ jointsatis ///
						terrdisp powerratio defense Rho
*Pull each variable, plus Rho
foreach i of num 7/14 {
	matrix m5_bl_sel[`i'-6,1] = beta[1,`i']
	matrix m5_bl_sel[`i'-6,2] = sqrt(cov[`i',`i'])
}
matrix m5_bl_sel[9,1] = beta[1,19]
matrix m5_bl_sel[9,2] = sqrt(cov[19,19])
***


* MID 3.0 War Model

heckprob war_mid3 joint_democ territory  jointsatis powerratio defense , ///
	sel(onset_mid3 = contiguous thom_rival min_min joint_democ  jointsatis ///
	terrdisp powerratio defense peaceyears pax2 pax3) cluster(dyadid)
eststo m6_bl
***
matrix m6_bl_out = J(5,2,.)
matrix coln m6_bl_out = coef se
matrix rown m6_bl_out = joint_democ territory jointsatis powerratio defense
matrix beta = e(b)
matrix cov = e(V)

foreach i of num 1/5 {
	matrix m6_bl_out[`i',1] = beta[1,`i']
	matrix m6_bl_out[`i',2] = sqrt(cov[`i',`i'])
}

matrix m6_bl_sel = J(9,2,.)
matrix coln m6_bl_sel = coef se
matrix rown m6_bl_sel = contiguous thom_rival min_min joint_democ jointsatis ///
						terrdisp powerratio defense Rho
*Pull each variable, plus Rho
foreach i of num 7/14 {
	matrix m6_bl_sel[`i'-6,1] = beta[1,`i']
	matrix m6_bl_sel[`i'-6,2] = sqrt(cov[`i',`i'])
}
matrix m6_bl_sel[9,1] = beta[1,19]
matrix m6_bl_sel[9,2] = sqrt(cov[19,19])
***


*************************************************************************************************************************************************

* ANALYSIS WITH GML 2.1 DATA


* We use version 2.1 of the GML data to be consistent with our analysis of Weeks (for which using version 2.1 was necessary to obtain incident-level data).
* Results using the GML V1.7 dataset that was used in the ISQ piece can be seen by substituting “gml1” for “gml2” in the regressions below.

* GML 2.1 Recip Model

heckprob recip_gml2 joint_democ territory jointsatis powerratio defense , ///
	sel(onset_gml2 = contiguous thom_rival min_min joint_democ  jointsatis ///
	terrdisp powerratio defense peaceyears pax2 pax3) cluster(dyadid)
eststo m1_bl
***
matrix m1_gml_out = J(5,2,.)
matrix coln m1_gml_out = coef se
matrix rown m1_gml_out = joint_democ territory jointsatis powerratio defense
matrix beta = e(b)
matrix cov = e(V)

foreach i of num 1/5 {
	matrix m1_gml_out[`i',1] = beta[1,`i']
	matrix m1_gml_out[`i',2] = sqrt(cov[`i',`i'])
}

matrix m1_gml_sel = J(9,2,.)
matrix coln m1_gml_sel = coef se
matrix rown m1_gml_sel = contiguous thom_rival min_min joint_democ jointsatis ///
						terrdisp powerratio defense Rho
*Pull each variable, plus Rho
foreach i of num 7/14 {
	matrix m1_gml_sel[`i'-6,1] = beta[1,`i']
	matrix m1_gml_sel[`i'-6,2] = sqrt(cov[`i',`i'])
}
matrix m1_gml_sel[9,1] = beta[1,19]
matrix m1_gml_sel[9,2] = sqrt(cov[19,19])
***

* GML 2.1 Use of Force Model

heckprob force_gml2 joint_democ territory  jointsatis powerratio defense , ///
	sel(onset_gml2 = contiguous thom_rival min_min joint_democ  jointsatis ///
	terrdisp powerratio defense peaceyears pax2 pax3) cluster(dyadid)
eststo m2_bl
***
matrix m2_gml_out = J(5,2,.)
matrix coln m2_gml_out = coef se
matrix rown m2_gml_out = joint_democ territory jointsatis powerratio defense
matrix beta = e(b)
matrix cov = e(V)

foreach i of num 1/5 {
	matrix m2_gml_out[`i',1] = beta[1,`i']
	matrix m2_gml_out[`i',2] = sqrt(cov[`i',`i'])
}

matrix m2_gml_sel = J(9,2,.)
matrix coln m2_gml_sel = coef se
matrix rown m2_gml_sel = contiguous thom_rival min_min joint_democ jointsatis ///
						terrdisp powerratio defense Rho
*Pull each variable, plus Rho
foreach i of num 7/14 {
	matrix m2_gml_sel[`i'-6,1] = beta[1,`i']
	matrix m2_gml_sel[`i'-6,2] = sqrt(cov[`i',`i'])
}
matrix m2_gml_sel[9,1] = beta[1,19]
matrix m2_gml_sel[9,2] = sqrt(cov[19,19])
***


* GML 2.1 Mutual Use of Force Model

heckprob mut_force_gml2 joint_democ territory  jointsatis powerratio defense  ///
	,sel(onset_gml2 = contiguous thom_rival min_min joint_democ  jointsatis ///
	terrdisp powerratio defense peaceyears pax2 pax3) cluster(dyadid)
eststo m3_bl
***
matrix m3_gml_out = J(5,2,.)
matrix coln m3_gml_out = coef se
matrix rown m3_gml_out = joint_democ territory jointsatis powerratio defense
matrix beta = e(b)
matrix cov = e(V)

foreach i of num 1/5 {
	matrix m3_gml_out[`i',1] = beta[1,`i']
	matrix m3_gml_out[`i',2] = sqrt(cov[`i',`i'])
}

matrix m3_gml_sel = J(9,2,.)
matrix coln m3_gml_sel = coef se
matrix rown m3_gml_sel = contiguous thom_rival min_min joint_democ jointsatis ///
						terrdisp powerratio defense Rho
*Pull each variable, plus Rho
foreach i of num 7/14 {
	matrix m3_gml_sel[`i'-6,1] = beta[1,`i']
	matrix m3_gml_sel[`i'-6,2] = sqrt(cov[`i',`i'])
}
matrix m3_gml_sel[9,1] = beta[1,19]
matrix m3_gml_sel[9,2] = sqrt(cov[19,19])
***


* GML 2.1 Fatal Model

heckprob fatal_gml2 joint_democ territory  jointsatis powerratio defense , ///
	sel(onset_gml2 = contiguous thom_rival min_min joint_democ  jointsatis ///
	terrdisp powerratio defense peaceyears pax2 pax3) cluster(dyadid)
eststo m4_bl
***
matrix m4_gml_out = J(5,2,.)
matrix coln m4_gml_out = coef se
matrix rown m4_gml_out = joint_democ territory jointsatis powerratio defense
matrix beta = e(b)
matrix cov = e(V)

foreach i of num 1/5 {
	matrix m4_gml_out[`i',1] = beta[1,`i']
	matrix m4_gml_out[`i',2] = sqrt(cov[`i',`i'])
}

matrix m4_gml_sel = J(9,2,.)
matrix coln m4_gml_sel = coef se
matrix rown m4_gml_sel = contiguous thom_rival min_min joint_democ jointsatis ///
						terrdisp powerratio defense Rho
*Pull each variable, plus Rho
foreach i of num 7/14 {
	matrix m4_gml_sel[`i'-6,1] = beta[1,`i']
	matrix m4_gml_sel[`i'-6,2] = sqrt(cov[`i',`i'])
}
matrix m4_gml_sel[9,1] = beta[1,19]
matrix m4_gml_sel[9,2] = sqrt(cov[19,19])
***


* GML 2.1 High Fatality Model

heckprob high_fatal_gml2 joint_democ territory  jointsatis powerratio defense , ///
	sel(onset_gml2 = contiguous thom_rival min_min joint_democ  jointsatis  ///
	terrdisp  powerratio defense peaceyears pax2 pax3) cluster(dyadid)
eststo m5_bl
***
matrix m5_gml_out = J(5,2,.)
matrix coln m5_gml_out = coef se
matrix rown m5_gml_out = joint_democ territory jointsatis powerratio defense
matrix beta = e(b)
matrix cov = e(V)

foreach i of num 1/5 {
	matrix m5_gml_out[`i',1] = beta[1,`i']
	matrix m5_gml_out[`i',2] = sqrt(cov[`i',`i'])
}

matrix m5_gml_sel = J(9,2,.)
matrix coln m5_gml_sel = coef se
matrix rown m5_gml_sel = contiguous thom_rival min_min joint_democ jointsatis ///
						terrdisp powerratio defense Rho
*Pull each variable, plus Rho
foreach i of num 7/14 {
	matrix m5_gml_sel[`i'-6,1] = beta[1,`i']
	matrix m5_gml_sel[`i'-6,2] = sqrt(cov[`i',`i'])
}
matrix m5_gml_sel[9,1] = beta[1,19]
matrix m5_gml_sel[9,2] = sqrt(cov[19,19])
***


* GML 2.1 War Model

heckprob war_gml2 joint_democ territory  jointsatis powerratio defense , ///
	sel(onset_gml2 = contiguous thom_rival min_min joint_democ  jointsatis ///
	terrdisp powerratio defense peaceyears pax2 pax3) cluster(dyadid)
eststo m6_bl
***
matrix m6_gml_out = J(5,2,.)
matrix coln m6_gml_out = coef se
matrix rown m6_gml_out = joint_democ territory jointsatis powerratio defense
matrix beta = e(b)
matrix cov = e(V)

foreach i of num 1/5 {
	matrix m6_gml_out[`i',1] = beta[1,`i']
	matrix m6_gml_out[`i',2] = sqrt(cov[`i',`i'])
}

matrix m6_gml_sel = J(9,2,.)
matrix coln m6_gml_sel = coef se
matrix rown m6_gml_sel = contiguous thom_rival min_min joint_democ jointsatis ///
						terrdisp powerratio defense Rho
*Pull each variable, plus Rho
foreach i of num 7/14 {
	matrix m6_gml_sel[`i'-6,1] = beta[1,`i']
	matrix m6_gml_sel[`i'-6,2] = sqrt(cov[`i',`i'])
}
matrix m6_gml_sel[9,1] = beta[1,19]
matrix m6_gml_sel[9,2] = sqrt(cov[19,19])
***



*************************************************************************************************************************************************

* ANALYSIS WITH MID 4.3 DATA

* MID 4.3 Recip Model

heckprob recip_mid4 joint_democ territory jointsatis powerratio defense , ///
	sel(onset_mid4 = contiguous thom_rival min_min joint_democ  jointsatis ///
	terrdisp powerratio defense peaceyears pax2 pax3) cluster(dyadid)
eststo m1_bl
***
matrix m1_mid4_out = J(5,2,.)
matrix coln m1_mid4_out = coef se
matrix rown m1_mid4_out = joint_democ territory jointsatis powerratio defense
matrix beta = e(b)
matrix cov = e(V)

foreach i of num 1/5 {
	matrix m1_mid4_out[`i',1] = beta[1,`i']
	matrix m1_mid4_out[`i',2] = sqrt(cov[`i',`i'])
}

matrix m1_mid4_sel = J(9,2,.)
matrix coln m1_mid4_sel = coef se
matrix rown m1_mid4_sel = contiguous thom_rival min_min joint_democ jointsatis ///
						terrdisp powerratio defense Rho
*Pull each variable, plus Rho
foreach i of num 7/14 {
	matrix m1_mid4_sel[`i'-6,1] = beta[1,`i']
	matrix m1_mid4_sel[`i'-6,2] = sqrt(cov[`i',`i'])
}
matrix m1_mid4_sel[9,1] = beta[1,19]
matrix m1_mid4_sel[9,2] = sqrt(cov[19,19])
***

* MID 4.3 Use of Force Model

heckprob force_mid4 joint_democ territory  jointsatis powerratio defense , ///
	sel(onset_mid4 = contiguous thom_rival min_min joint_democ  jointsatis ///
	terrdisp powerratio defense peaceyears pax2 pax3) cluster(dyadid)
eststo m2_bl
***
matrix m2_mid4_out = J(5,2,.)
matrix coln m2_mid4_out = coef se
matrix rown m2_mid4_out = joint_democ territory jointsatis powerratio defense
matrix beta = e(b)
matrix cov = e(V)

foreach i of num 1/5 {
	matrix m2_mid4_out[`i',1] = beta[1,`i']
	matrix m2_mid4_out[`i',2] = sqrt(cov[`i',`i'])
}

matrix m2_mid4_sel = J(9,2,.)
matrix coln m2_mid4_sel = coef se
matrix rown m2_mid4_sel = contiguous thom_rival min_min joint_democ jointsatis ///
						terrdisp powerratio defense Rho
*Pull each variable, plus Rho
foreach i of num 7/14 {
	matrix m2_mid4_sel[`i'-6,1] = beta[1,`i']
	matrix m2_mid4_sel[`i'-6,2] = sqrt(cov[`i',`i'])
}
matrix m2_mid4_sel[9,1] = beta[1,19]
matrix m2_mid4_sel[9,2] = sqrt(cov[19,19])
***


* MID 4.3 Mutual Use of Force Model

heckprob mut_force_mid4 joint_democ territory  jointsatis powerratio defense  ///
	,sel(onset_mid4 = contiguous thom_rival min_min joint_democ  jointsatis ///
	terrdisp powerratio defense peaceyears pax2 pax3) cluster(dyadid)
eststo m3_bl
***
matrix m3_mid4_out = J(5,2,.)
matrix coln m3_mid4_out = coef se
matrix rown m3_mid4_out = joint_democ territory jointsatis powerratio defense
matrix beta = e(b)
matrix cov = e(V)

foreach i of num 1/5 {
	matrix m3_mid4_out[`i',1] = beta[1,`i']
	matrix m3_mid4_out[`i',2] = sqrt(cov[`i',`i'])
}

matrix m3_mid4_sel = J(9,2,.)
matrix coln m3_mid4_sel = coef se
matrix rown m3_mid4_sel = contiguous thom_rival min_min joint_democ jointsatis ///
						terrdisp powerratio defense Rho
*Pull each variable, plus Rho
foreach i of num 7/14 {
	matrix m3_mid4_sel[`i'-6,1] = beta[1,`i']
	matrix m3_mid4_sel[`i'-6,2] = sqrt(cov[`i',`i'])
}
matrix m3_mid4_sel[9,1] = beta[1,19]
matrix m3_mid4_sel[9,2] = sqrt(cov[19,19])
***


* MID 4.3 Fatal Model

heckprob fatal_mid4 joint_democ territory  jointsatis powerratio defense , ///
	sel(onset_mid4 = contiguous thom_rival min_min joint_democ  jointsatis ///
	terrdisp powerratio defense peaceyears pax2 pax3) cluster(dyadid)
eststo m4_bl
***
matrix m4_mid4_out = J(5,2,.)
matrix coln m4_mid4_out = coef se
matrix rown m4_mid4_out = joint_democ territory jointsatis powerratio defense
matrix beta = e(b)
matrix cov = e(V)

foreach i of num 1/5 {
	matrix m4_mid4_out[`i',1] = beta[1,`i']
	matrix m4_mid4_out[`i',2] = sqrt(cov[`i',`i'])
}

matrix m4_mid4_sel = J(9,2,.)
matrix coln m4_mid4_sel = coef se
matrix rown m4_mid4_sel = contiguous thom_rival min_min joint_democ jointsatis ///
						terrdisp powerratio defense Rho
*Pull each variable, plus Rho
foreach i of num 7/14 {
	matrix m4_mid4_sel[`i'-6,1] = beta[1,`i']
	matrix m4_mid4_sel[`i'-6,2] = sqrt(cov[`i',`i'])
}
matrix m4_mid4_sel[9,1] = beta[1,19]
matrix m4_mid4_sel[9,2] = sqrt(cov[19,19])
***


* MID 4.3 High Fatality Model

heckprob high_fatal_mid4 joint_democ territory  jointsatis powerratio defense , ///
	sel(onset_mid4 = contiguous thom_rival min_min joint_democ  jointsatis  ///
	terrdisp  powerratio defense peaceyears pax2 pax3) cluster(dyadid)
eststo m5_bl
***
matrix m5_mid4_out = J(5,2,.)
matrix coln m5_mid4_out = coef se
matrix rown m5_mid4_out = joint_democ territory jointsatis powerratio defense
matrix beta = e(b)
matrix cov = e(V)

foreach i of num 1/5 {
	matrix m5_mid4_out[`i',1] = beta[1,`i']
	matrix m5_mid4_out[`i',2] = sqrt(cov[`i',`i'])
}

matrix m5_mid4_sel = J(9,2,.)
matrix coln m5_mid4_sel = coef se
matrix rown m5_mid4_sel = contiguous thom_rival min_min joint_democ jointsatis ///
						terrdisp powerratio defense Rho
*Pull each variable, plus Rho
foreach i of num 7/14 {
	matrix m5_mid4_sel[`i'-6,1] = beta[1,`i']
	matrix m5_mid4_sel[`i'-6,2] = sqrt(cov[`i',`i'])
}
matrix m5_mid4_sel[9,1] = beta[1,19]
matrix m5_mid4_sel[9,2] = sqrt(cov[19,19])
***


* MID 4.3 War Model

heckprob war_mid4 joint_democ territory  jointsatis powerratio defense , ///
	sel(onset_mid4 = contiguous thom_rival min_min joint_democ  jointsatis ///
	terrdisp powerratio defense peaceyears pax2 pax3) cluster(dyadid)
eststo m6_bl
***
matrix m6_mid4_out = J(5,2,.)
matrix coln m6_mid4_out = coef se
matrix rown m6_mid4_out = joint_democ territory jointsatis powerratio defense
matrix beta = e(b)
matrix cov = e(V)

foreach i of num 1/5 {
	matrix m6_mid4_out[`i',1] = beta[1,`i']
	matrix m6_mid4_out[`i',2] = sqrt(cov[`i',`i'])
}

matrix m6_mid4_sel = J(9,2,.)
matrix coln m6_mid4_sel = coef se
matrix rown m6_mid4_sel = contiguous thom_rival min_min joint_democ jointsatis ///
						terrdisp powerratio defense Rho
*Pull each variable, plus Rho
foreach i of num 7/14 {
	matrix m6_mid4_sel[`i'-6,1] = beta[1,`i']
	matrix m6_mid4_sel[`i'-6,2] = sqrt(cov[`i',`i'])
}
matrix m6_mid4_sel[9,1] = beta[1,19]
matrix m6_mid4_sel[9,2] = sqrt(cov[19,19])
***

*******************************************************************************************************************************

* MAKE FIGURES

coefplot ///
	(matrix(m1_bl_sel[.,1]), se(m1_mid4_sel[.,2]) m(smcircle) mlw(medium) msize(medlarge) mfcolor("253 174 97*.4") mlc("253 174 97")  ciop(lcolor("253 174 97" "253 174 97") lwidth(medthick thick) ) label(MID 3.0)) ///
	(matrix(m1_gml_sel[.,1]), se(m1_mid4_sel[.,2])  m(smcircle) mlw(medium) msize(medlarge) mfcolor("215 25 28*.4")  mlc("215 25 28") ciop(lcolor("215 25 28" "215 25 28") lwidth(medthick thick) ) label(GML Data)) ///
	(matrix(m1_mid4_sel[.,1]), se(m1_mid4_sel[.,2])  m(smcircle) mlw(medium) msize(medlarge) mfcolor("31 120 180*.4")  mlc("31 120 180") ciop(lcolor("31 120 180" "31 120 180") lwidth(medthick thick) ) label(MID 4.3)) ///
	,  scheme(s1mono)  ///
	coeflabels()   ///
	legend(col(3) ring()  region(lcolor(none)) size(2)) ///
	xtitle("Coefficient Value" , size(small)) ///
	ytitle(" ", size(small) margin(b=0 l=0 t=0 r=0) ) ///
	subtitle("Onset" , size(small)) ///
	ysca(noline titlegap(0)) ///
	xsca( noextend titlegap(3) range(-1.25,1.25)) ///
	xlabel(-1 -.5 0 .5 1 , labsize(small) angle(0)  ) ///
	ylabel(,notick   labsize(small)  ) ///
	plotregion(color(gs15) style(none)) ///
	xline(0, lcolor(gray) lpattern(dash)) ///
	ysize(20) xsize() fysize() ///
		graphregion(margin(l=0 r=-2)) ///
	levels(95 90) saving(m1_sel, replace) 
coefplot ///
	(matrix(m1_bl_out[.,1]), se(m1_bl_out[.,2]) m(smcircle) mlw(medium) msize(medium) mfcolor("253 174 97*.4") mlc("253 174 97")  ciop(lcolor("253 174 97" "253 174 97") lwidth(medthick thick) ) label(MID 3.0)) ///
	(matrix(m1_gml_out[.,1]), se(m1_gml_out[.,2])  m(smcircle) mlw(medium) msize(medium) mfcolor("215 25 28*.4")  mlc("215 25 28") ciop(lcolor("215 25 28" "215 25 28") lwidth(medthick thick) ) label(GML Data)) ///
	(matrix(m1_mid4_out[.,1]), se(m1_mid4_out[.,2])  m(smcircle) mlw(medium) msize(medium) mfcolor("31 120 180*.4")  mlc("31 120 180") ciop(lcolor("31 120 180" "31 120 180") lwidth(medthick thick) ) label(MID 4.3)) ///
	,  scheme(s1mono)  ///
	coeflabels()   ///
	legend(off) ///
	xtitle("") ///
	ytitle(" ", size(small) margin(b=0 l=5.5 t=0 r=0)) ///
	title("Model 1"  , size(medsmall)) ///
	subtitle("Reciprocation" , size(small)) ///
	ysca(noline titlegap(0)) ///
	xsca( noline noextend titlegap(3) range(-1.25,1.25)) ///
	xlabel(none  ) ///
	ylabel(,notick   labsize(small)  ) ///
	plotregion(color(gs15) style(none)) ///
	xline(0, lcolor(gray) lpattern(dash)) ///
	ysize(5) xsize() fysize(80) ///
	graphregion(margin(l=0 r=-2 b=-8)) ///
	levels(95 90) saving(m1_out, replace)
graph combine m1_out.gph m1_sel.gph,  xsize(3) ysize(6) iscale(1.2) ///
	 scheme(s1mono) rows(2) saving(m1, replace) 

	 
	 
coefplot ///
	(matrix(m2_bl_sel[.,1]), se(m2_mid4_sel[.,2]) m(smcircle) mlw(medium) msize(medlarge) mfcolor("253 174 97*.4") mlc("253 174 97")  ciop(lcolor("253 174 97" "253 174 97") lwidth(medthick thick) ) label(MID 3.0)) ///
	(matrix(m2_gml_sel[.,1]), se(m2_mid4_sel[.,2])  m(smcircle) mlw(medium) msize(medlarge) mfcolor("215 25 28*.4")  mlc("215 25 28") ciop(lcolor("215 25 28" "215 25 28") lwidth(medthick thick) ) label(GML Data)) ///
	(matrix(m2_mid4_sel[.,1]), se(m2_mid4_sel[.,2])  m(smcircle) mlw(medium) msize(medlarge) mfcolor("31 120 180*.4")  mlc("31 120 180") ciop(lcolor("31 120 180" "31 120 180") lwidth(medthick thick) ) label(MID 4.3)) ///
	,  scheme(s1mono)  ///
	coeflabels()   ///
	legend(col(3) ring()  region(lcolor(none)) size(2)) ///
	xtitle("Coefficient Value" , size(small)) ///
	ytitle(" ", size(small) margin(0 15 0 0) ) ///
	subtitle("Onset" , size(small)) ///
	ysca(noline titlegap(0)) ///
	xsca( noextend titlegap(3) range(-1.25,1.25)) ///
	xlabel(-1 -.5 0 .5 1 , labsize(small) angle(0)  ) ///
	ylabel(none ) ///
	plotregion(color(gs15) style(none)) ///
	xline(0, lcolor(gray) lpattern(dash)) ///
	ysize(20) xsize() fysize() ///
	graphregion(margin(l=0 r=18)) ///
	levels(95 90) saving(m2_sel, replace) 
coefplot ///
	(matrix(m2_bl_out[.,1]), se(m2_bl_out[.,2]) m(smcircle) mlw(medium) msize(medium) mfcolor("253 174 97*.4") mlc("253 174 97")  ciop(lcolor("253 174 97" "253 174 97") lwidth(medthick thick) ) label(MID 3.0)) ///
	(matrix(m2_gml_out[.,1]), se(m2_gml_out[.,2])  m(smcircle) mlw(medium) msize(medium) mfcolor("215 25 28*.4")  mlc("215 25 28") ciop(lcolor("215 25 28" "215 25 28") lwidth(medthick thick) ) label(GML Data)) ///
	(matrix(m2_mid4_out[.,1]), se(m2_mid4_out[.,2])  m(smcircle) mlw(medium) msize(medium) mfcolor("31 120 180*.4")  mlc("31 120 180") ciop(lcolor("31 120 180" "31 120 180") lwidth(medthick thick) ) label(MID 4.3)) ///
	,  scheme(s1mono)  ///
	coeflabels()   ///
	legend(off) ///
	xtitle("") ///
	ytitle(" ", size(small) margin(0 15 0 0)) ///
	title("Model 2"  , size(medsmall)) ///
	subtitle("Use of Force" , size(small)) ///
	ysca(noline titlegap(0)) ///
	xsca( noline noextend titlegap(3) range(-1.25,1.25)) ///
	xlabel(none) ///
	ylabel(none) ///
	plotregion(color(gs15) style(none)) ///
	xline(0, lcolor(gray) lpattern(dash)) ///
	ysize(5) xsize() fysize(80) ///
	graphregion(margin(l=0 r=18 b=-8)) ///
	levels(95 90) saving(m2_out, replace)

graph combine m2_out.gph m2_sel.gph,  xsize(3) ysize(6) iscale(1.2) ///
	 scheme(s1mono) rows(2) saving(m2, replace) 
	 
	 
*Model 3
coefplot ///
	(matrix(m3_bl_sel[.,1]), se(m3_mid4_sel[.,2]) m(smcircle) mlw(medium) msize(medlarge) mfcolor("253 174 97*.4") mlc("253 174 97")  ciop(lcolor("253 174 97" "253 174 97") lwidth(medthick thick) ) label(MID 3.0)) ///
	(matrix(m3_gml_sel[.,1]), se(m3_mid4_sel[.,2])  m(smcircle) mlw(medium) msize(medlarge) mfcolor("215 25 28*.4")  mlc("215 25 28") ciop(lcolor("215 25 28" "215 25 28") lwidth(medthick thick) ) label(GML Data)) ///
	(matrix(m3_mid4_sel[.,1]), se(m3_mid4_sel[.,2])  m(smcircle) mlw(medium) msize(medlarge) mfcolor("31 120 180*.4")  mlc("31 120 180") ciop(lcolor("31 120 180" "31 120 180") lwidth(medthick thick) ) label(MID 4.3)) ///
	,  scheme(s1mono)  ///
	coeflabels()   ///
	legend(col(3) ring()  region(lcolor(none)) size(2)) ///
	xtitle("Coefficient Value" , size(small)) ///
	ytitle(" ", size(small) ) ///
	subtitle("Onset" , size(small)) ///
	ysca(noline titlegap(0)) ///
	xsca( noextend titlegap(3) range(-1.25,1.25)) ///
	xlabel(-1 -.5 0 .5 1 , labsize(small) angle(0)  ) ///
	ylabel(none ) ///
	plotregion(color(gs15) style(none)) ///
	xline(0, lcolor(gray) lpattern(dash)) ///
	ysize(20) xsize() fysize() ///
	graphregion(margin(l=9 r=23.5)) ///
	levels(95 90) saving(m3_sel, replace) 
coefplot ///
	(matrix(m3_bl_out[.,1]), se(m3_bl_out[.,2]) m(smcircle) mlw(medium) msize(medium) mfcolor("253 174 97*.4") mlc("253 174 97")  ciop(lcolor("253 174 97" "253 174 97") lwidth(medthick thick) ) label(MID 3.0)) ///
	(matrix(m3_gml_out[.,1]), se(m3_gml_out[.,2])  m(smcircle) mlw(medium) msize(medium) mfcolor("215 25 28*.4")  mlc("215 25 28") ciop(lcolor("215 25 28" "215 25 28") lwidth(medthick thick) ) label(GML Data)) ///
	(matrix(m3_mid4_out[.,1]), se(m3_mid4_out[.,2])  m(smcircle) mlw(medium) msize(medium) mfcolor("31 120 180*.4")  mlc("31 120 180") ciop(lcolor("31 120 180" "31 120 180") lwidth(medthick thick) ) label(MID 4.3)) ///
	,  scheme(s1mono)  ///
	coeflabels()   ///
	legend(off) ///
	xtitle("") ///
	ytitle(" ", size(small) margin(0 0 0 0)) ///
	title("Model 3"  , size(medsmall)) ///
	subtitle("Mutual Use of Force" , size(small)) ///
	ysca(noline titlegap(0)) ///
	xsca( noline noextend titlegap(3) range(-1.25,1.25)) ///
	xlabel(none  ) ///
	ylabel(none ) ///
	plotregion(color(gs15) style(none)) ///
	xline(0, lcolor(gray) lpattern(dash)) ///
	ysize(5) xsize() fysize(80) ///
	graphregion(margin(l=9 r=23.5 b=-8)) ///
	levels(95 90) saving(m3_out, replace)
	 
graph combine m3_out.gph m3_sel.gph,  xsize(3) ysize(6) iscale(1.2) ///
	 scheme(s1mono) rows(2) saving(m3, replace) 
	

*Model4
coefplot ///
	(matrix(m4_bl_sel[.,1]), se(m4_mid4_sel[.,2]) m(smcircle) mlw(medium) msize(medlarge) mfcolor("253 174 97*.4") mlc("253 174 97")  ciop(lcolor("253 174 97" "253 174 97") lwidth(medthick thick) ) label(MID 3.0)) ///
	(matrix(m4_gml_sel[.,1]), se(m4_mid4_sel[.,2])  m(smcircle) mlw(medium) msize(medlarge) mfcolor("215 25 28*.4")  mlc("215 25 28") ciop(lcolor("215 25 28" "215 25 28") lwidth(medthick thick) ) label(GML Data)) ///
	(matrix(m4_mid4_sel[.,1]), se(m4_mid4_sel[.,2])  m(smcircle) mlw(medium) msize(medlarge) mfcolor("31 120 180*.4")  mlc("31 120 180") ciop(lcolor("31 120 180" "31 120 180") lwidth(medthick thick) ) label(MID 4.3)) ///
	,  scheme(s1mono)  ///
	coeflabels()   ///
	legend(col(3) ring()  region(lcolor(none)) size(2)) ///
	xtitle("Coefficient Value" , size(small)) ///
	ytitle(" ", size(small) ) ///
	subtitle("Onset" , size(small)) ///
	ysca(noline titlegap(0)) ///
	xsca( noextend titlegap(3) range(-1.25,1.25)) ///
	xlabel(-1 -.5 0 .5 1 , labsize(small) angle(0)  ) ///
	ylabel(,notick   labsize(small)  ) ///
	plotregion(color(gs15) style(none)) ///
	xline(0, lcolor(gray) lpattern(dash)) ///
	ysize(20) xsize() fysize() ///
	graphregion(margin(l=0 r=-2)) ///
	levels(95 90) saving(m4_sel, replace) 
coefplot ///
	(matrix(m4_bl_out[.,1]), se(m4_bl_out[.,2]) m(smcircle) mlw(medium) msize(medium) mfcolor("253 174 97*.4") mlc("253 174 97")  ciop(lcolor("253 174 97" "253 174 97") lwidth(medthick thick) ) label(MID 3.0)) ///
	(matrix(m4_gml_out[.,1]), se(m4_gml_out[.,2])  m(smcircle) mlw(medium) msize(medium) mfcolor("215 25 28*.4")  mlc("215 25 28") ciop(lcolor("215 25 28" "215 25 28") lwidth(medthick thick) ) label(GML Data)) ///
	(matrix(m4_mid4_out[.,1]), se(m4_mid4_out[.,2])  m(smcircle) mlw(medium) msize(medium) mfcolor("31 120 180*.4")  mlc("31 120 180") ciop(lcolor("31 120 180" "31 120 180") lwidth(medthick thick) ) label(MID 4.3)) ///
	,  scheme(s1mono)  ///
	coeflabels()   ///
	legend(off) ///
	xtitle("") ///
	ytitle(" ", size(small) margin(0 5.5 0 0)) ///
	title("Model 4"  , size(medsmall)) ///
	subtitle("Fatal MID" , size(small)) ///
	ysca(noline titlegap(0)) ///
	xsca( noline noextend titlegap(3) range(-1.25,1.25)) ///
	xlabel(none  ) ///
	ylabel(,notick   labsize(small)  ) ///
	plotregion(color(gs15) style(none)) ///
	xline(0, lcolor(gray) lpattern(dash)) ///
	ysize(5) xsize() fysize(80) ///
	graphregion(margin(l=0 r=-2 b=-8)) ///
	levels(95 90) saving(m4_out, replace)
graph combine m4_out.gph m4_sel.gph,  xsize(3) ysize(6) iscale(1.2) ///
	 scheme(s1mono) rows(2) saving(m4, replace) 
	
	
	
*Model 5
coefplot ///
	(matrix(m5_bl_sel[.,1]), se(m5_mid4_sel[.,2]) m(smcircle) mlw(medium) msize(medlarge) mfcolor("253 174 97*.4") mlc("253 174 97")  ciop(lcolor("253 174 97" "253 174 97") lwidth(medthick thick) ) label(MID 3.0)) ///
	(matrix(m5_gml_sel[.,1]), se(m5_mid4_sel[.,2])  m(smcircle) mlw(medium) msize(medlarge) mfcolor("215 25 28*.4")  mlc("215 25 28") ciop(lcolor("215 25 28" "215 25 28") lwidth(medthick thick) ) label(GML Data)) ///
	(matrix(m5_mid4_sel[.,1]), se(m5_mid4_sel[.,2])  m(smcircle) mlw(medium) msize(medlarge) mfcolor("31 120 180*.4")  mlc("31 120 180") ciop(lcolor("31 120 180" "31 120 180") lwidth(medthick thick) ) label(MID 4.3)) ///
	,  scheme(s1mono)  ///
	coeflabels()   ///
	legend(col(3) ring()  region(lcolor(none)) size(2)) ///
	xtitle("Coefficient Value" , size(small)) ///
	ytitle(" ", size(small) margin(0 15 0 0) ) ///
	subtitle("Onset" , size(small)) ///
	ysca(noline titlegap(0)) ///
	xsca( noextend titlegap(3) range(-1.25,1.25)) ///
	xlabel(-1 -.5 0 .5 1 , labsize(small) angle(0)  ) ///
	ylabel(none ) ///
	plotregion(color(gs15) style(none)) ///
	xline(0, lcolor(gray) lpattern(dash)) ///
	ysize(20) xsize() fysize() ///
	graphregion(margin(l=0 r=18)) ///
	levels(95 90) saving(m5_sel, replace) 
coefplot ///
	(matrix(m5_bl_out[.,1]), se(m5_bl_out[.,2]) m(smcircle) mlw(medium) msize(medium) mfcolor("253 174 97*.4") mlc("253 174 97")  ciop(lcolor("253 174 97" "253 174 97") lwidth(medthick thick) ) label(MID 3.0)) ///
	(matrix(m5_gml_out[.,1]), se(m5_gml_out[.,2])  m(smcircle) mlw(medium) msize(medium) mfcolor("215 25 28*.4")  mlc("215 25 28") ciop(lcolor("215 25 28" "215 25 28") lwidth(medthick thick) ) label(GML Data)) ///
	(matrix(m5_mid4_out[.,1]), se(m5_mid4_out[.,2])  m(smcircle) mlw(medium) msize(medium) mfcolor("31 120 180*.4")  mlc("31 120 180") ciop(lcolor("31 120 180" "31 120 180") lwidth(medthick thick) ) label(MID 4.3)) ///
	,  scheme(s1mono)  ///
	coeflabels()   ///
	legend(off) ///
	xtitle("") ///
	ytitle(" ", size(small) margin(0 15 0 0)) ///
	title("Model 5"  , size(medsmall)) ///
	subtitle("MID with over 250 Fatalities" , size(small)) ///
	ysca(noline titlegap(0)) ///
	xsca( noline noextend titlegap(3) range(-1.25,1.25)) ///
	xlabel(none  ) ///
	ylabel(none) ///
	plotregion(color(gs15) style(none)) ///
	xline(0, lcolor(gray) lpattern(dash)) ///
	ysize(5) xsize() fysize(80) ///
	graphregion(margin(l=0 r=18 b=-8)) ///
	levels(95 90) saving(m5_out, replace)
graph combine m5_out.gph m5_sel.gph,  xsize(3) ysize(6) iscale(1.2) ///
	 scheme(s1mono) rows(2) saving(m5, replace) 

*Model 6
coefplot ///
	(matrix(m6_bl_sel[.,1]), se(m6_mid4_sel[.,2]) m(smcircle) mlw(medium) msize(medlarge) mfcolor("253 174 97*.4") mlc("253 174 97")  ciop(lcolor("253 174 97" "253 174 97") lwidth(medthick thick) ) label(MID 3.0)) ///
	(matrix(m6_gml_sel[.,1]), se(m6_mid4_sel[.,2])  m(smcircle) mlw(medium) msize(medlarge) mfcolor("215 25 28*.4")  mlc("215 25 28") ciop(lcolor("215 25 28" "215 25 28") lwidth(medthick thick) ) label(GML Data)) ///
	(matrix(m6_mid4_sel[.,1]), se(m6_mid4_sel[.,2])  m(smcircle) mlw(medium) msize(medlarge) mfcolor("31 120 180*.4")  mlc("31 120 180") ciop(lcolor("31 120 180" "31 120 180") lwidth(medthick thick) ) label(MID 4.3)) ///
	,  scheme(s1mono)  ///
	coeflabels()   ///
	legend(col(3) ring()  region(lcolor(none)) size(2)) ///
	xtitle("Coefficient Value" , size(small)) ///
	ytitle(" ", size(small) ) ///
	subtitle("Onset" , size(small)) ///
	ysca(noline titlegap(0)) ///
	xsca( noextend titlegap(3) range(-1,2)) ///
	xlabel(-1 -.5 0 .5 1 1.5 2 , labsize(small) angle(0)  ) ///
	ylabel(none ) ///
	plotregion(color(gs15) style(none)) ///
	xline(0, lcolor(gray) lpattern(dash)) ///
	ysize(20) xsize() fysize() ///
	graphregion(margin(l=9 r=23.5)) ///
	levels(95 90) saving(m6_sel, replace) 
coefplot ///
	(matrix(m6_bl_out[.,1]), se(m6_bl_out[.,2]) m(smcircle) mlw(medium) msize(medium) mfcolor("253 174 97*.4") mlc("253 174 97")  ciop(lcolor("253 174 97" "253 174 97") lwidth(medthick thick) ) label(MID 3.0)) ///
	(matrix(m6_gml_out[.,1]), se(m6_gml_out[.,2])  m(smcircle) mlw(medium) msize(medium) mfcolor("215 25 28*.4")  mlc("215 25 28") ciop(lcolor("215 25 28" "215 25 28") lwidth(medthick thick) ) label(GML Data)) ///
	(matrix(m6_mid4_out[.,1]), se(m6_mid4_out[.,2])  m(smcircle) mlw(medium) msize(medium) mfcolor("31 120 180*.4")  mlc("31 120 180") ciop(lcolor("31 120 180" "31 120 180") lwidth(medthick thick) ) label(MID 4.3)) ///
	,  scheme(s1mono)  ///
	coeflabels()   ///
	legend(off) ///
	xtitle("") ///
	ytitle(" ", size(small) margin(0 0 0 0)) ///
	title("Model 6"  , size(medsmall)) ///
	subtitle("War" , size(small)) ///
	ysca(noline titlegap(0)) ///
	xsca( noline noextend titlegap(3) range(-1,2)) ///
	xlabel(none  ) ///
	ylabel(none ) ///
	plotregion(color(gs15) style(none)) ///
	xline(0, lcolor(gray) lpattern(dash)) ///
	ysize(5) xsize() fysize(80) ///
	graphregion(margin(l=9 r=23.5 b=-8)) ///
	levels(95 90) saving(m6_out, replace)
graph combine m6_out.gph m6_sel.gph,  xsize(3) ysize(6) iscale(1.2) ///
	 scheme(s1mono) rows(2) saving(m6, replace) 
	
